Bioenergetic control of soil carbon dynamics across depth

Soil carbon dynamics is strongly controlled by depth globally, with increasingly slow dynamics found at depth. The mechanistic basis remains however controversial, limiting our ability to predict carbon cycle-climate feedbacks. Here we combine radiocarbon and thermal analyses with long-term incubations in absence/presence of continuously 13C/14C-labelled plants to show that bioenergetic constraints of decomposers consistently drive the depth-dependency of soil carbon dynamics over a range of mineral reactivity contexts. The slow dynamics of subsoil carbon is tightly related to both its low energy density and high activation energy of decomposition, leading to an unfavourable ‘return-on-energy-investment’ for decomposers. We also observe strong acceleration of millennia-old subsoil carbon decomposition induced by roots (‘rhizosphere priming’), showing that sufficient supply of energy by roots is able to alleviate the strong energy limitation of decomposition. These findings demonstrate that subsoil carbon persistence results from its poor energy quality together with the lack of energy supply by roots due to their low density at depth.


Soil biogeochemical properties Radiocarbon ( 14 C) measurements and SOC turnover time estimation
Once in the 14 C laboratory, soil samples were crushed under 200 µm to homogeneize it. As no carbonate was expected, samples were not acidified prior to the measurement. According to SOC content of the samples, 14 C measurements were performed using either a solid source for topsoil samples or gas source for subsoil samples. Aliquotes were sampled in tin capsules to get either 1 mg of SOC for topsoils or twice 70 µg of SOC for subsoils. For topsoil samples, CO2 evolved from SOC was graphitized using an automated AGE 3 graphitization device 1 , and 14 C measurement was performed on ECHoMICADAS 2-4 through the solid source. For each subsoil samples, two 14 C measurements were performed on ECHoMICADAS through the gaz source connected by a Gas Ion Source (GIS) 5 to an elementar analyser (EA) evolving SOC to CO2. The 14 C signature of root biomass were performed as decribed above for topsoil samples. Barium carbonate from soil respired CO2 was transformed into CO2 in a semi-automated carbonate line 6 . Evolved CO2 was sealed under vacuum in pyrex microtube and introduced into the ECHoMICADAS gas source through the cracking-GIS interface 5 . All 14 C results have been corrected for mass-dependent isotopic fractionation using AMS-derived 13 C measurements, and were expressed as deviations from the absolute (decay-corrected) Oxalic Acid I (OX1) standard (Δ 14 C, in ‰) 7,8 . The average measurement precision of the Δ 14 C values was 2.4 ‰.
We estimated SOC turnover time based on radiocarbon measurements using a modelling approach. The following time-dependent, homogeneous one-pool model 7,9 was used: given 14 where at time t, F 14C SOC is the 14 C content of SOC, SOC is the SOC stock, I is the rate of C input from the atmosphere to SOC, F 14C atm is the 14 C content of CO2 in the local atmosphere, TR is the mean C transit time through living plant material, τ is the mean SOC turnover time and λ is the radioactive decay constant for 14 C (1.21 × 10 −4 year −1 ). Assuming TR = 1 year and SOC stocks to be at steady-state so that Ct = Ct-1 = I × τ, the equation (1) reduces to: (1 − 1 − ) (2) where t-1 is the year preceding time t. We ran the model from 50 kyr BP until the year 2016 using the 'SoilR' package to calculate the predicted SOC Δ 14 C at the year of sampling for a range of τ values (1 to 30,000 years) 10 . We used Δ 14 C atmospheric values from Reimer et al., (2020) 11 for the period 0-50 000 years BP, from Hua et al., (2013) 12 for the period 1950-2010, and we calculated extrapolated values using an exponential smoothing state-space modeling approach for the period 2010-2016 10 . We then derived τ values from our Δ 14 C measurements for each sample based on the relationship between τ and predicted Δ 14 C (Supplementary Fig.  6).

Thermal analyses
Rock-Eval® thermal analysis consisting in evolved gas analysis during ramped combustion was performed using a Rock-Eval® 6 Turbo device (Vinci Technologies, France) following a procedure adapted for SOM analysis 13 to measure the activation energy (Ea) of thermal SOC decomposition. Briefly, ca. 60 mg of ground (<250 μm) samples were subjected to sequential pyrolysis and oxidation phases. The pyrolysis phase was carried out in an N2 atmosphere with a 3 min isotherm at 200 °C followed by a temperature ramp from 200 to 650 °C at a heating rate of 30 °C min -1 . The oxidation phase was carried out in laboratory air atmosphere with a 1 min isotherm at 300 °C followed by a temperature ramp from 300 to 850 °C at a heating rate of 20 °C min -1 and a final 5 min isotherm at 850 °C. Hydrocarbon effluents (HxCy) were quantified by flame ionization detection during the pyrolysis phase, while CO and CO2 were quantified by infrared detection during both ramping phases. Each Rock-Eval® thermal analysis generated five thermograms corresponding to hydrocarbon effluents (HxCy-pyrolysis thermogram), CO (CO-pyrolysis thermogram) and CO2 (CO2-pyrolysis thermogram) measured at each second during the pyrolysis phase, and to the CO (CO-oxidation thermogram) and CO2 (CO2-oxidation thermogram) measured at each second during the oxidation phase. Thermograms were integrated on different time intervals depending on the thermogram. The integration omitted the first 200 seconds of the analysis for the three thermograms of the pyrolysis phase. The integration ended at the time of analysis corresponding to the maximum oven temperatures of 650 °C (HxCy-pyrolysis thermogram), 560 °C (CO-pyrolysis and CO2-pyrolysis thermograms), 850 °C (CO-oxidation thermogram) and 611 °C (CO2-oxidation thermogram). These intervals of integration prevented any interference by inorganic carbon from most soil carbonates 14 . Before determination of the activation energy (Ea) of SOC decomposition, the three thermograms of the pyrolysis phase and the two thermograms of the oxidation phase were combined into single thermograms of mass-equivalent C evolved during each ramping phase 15 .
Assuming first-order reaction kinetics during ramped combustion, we used a regularized, inverse method to determine the continuous distribution of Ea that best predicts the measured SOC decay profile 16 . Thermograms of mass-equivalent C evolved during both pyrolysis and oxidation were analysed separately using the 'rampedpyrox' Python package 17 , resulting in distinct Ea distributions for each ramping phase. For both the thermograms and Ea distributions, the data from each phase were then merged by averaging weighted based on the relative amount of C evolved during each ramping phase. Each combined continuous Ea distribution was then integrated to calculate the mean (µEa) and standard deviation (σEa) of activation energy (kJ mol -1 SOM) 16 . The Ea represents the energy input required for SOC combustion and was used here as a proxy for the energetic barriers to microbial SOC decomposition, that is the energy input needed from microbes producing exoenzymes to access and metabolize it 18 . We acknowledge that these computed Ea values of thermal SOC decomposition are expected to be much higher than for naturally occurring SOC biodegradation catalysed by microbial exoenzymes 19 . However, many studies combining radiocarbon and thermal analyses have found Ea of SOC combustion to be strongly related to SOC biogeochemical stability, with increasing Ea for SOC having increasing SOC radiocarbon ages [20][21][22] . We thus assumed here that the Ea of SOC combustion remains a good proxy of the energy investment in exoenzymes production needed for microbial decomposers to acquire this SOC. Though the pyrolysis phase may have overestimated our Ea values due to charring effects 23 , it has been found that these charring reactions do not affect the determination and interpretations of relative thermal stability of SOC 24 . Previous studies further showed that ramped pyrolysis and oxidation yielded similar Ea results 20,21 .
Two standard Rock-Eval® 6 parameters describing SOM bulk chemistry in term of hydrogen and oxygen composition were determined 14 , that are the hydrogen and oxygen indices (HI and OI). The HI index was calculated the amount of hydrocarbons (HxCy) formed during thermal pyrolysis of the sample between 200 and 650 °C divided by the total SOC of the sample. The OI was calculated using the following equation 25 27 . A preliminary study measuring the elemental H:C and O:C ratios of SOM based on nuclear magnetic resonance spectroscopy coupled to a molecular mixing model found strong relations with respectively HI (H:C = 1.21 + HI × 8.20 × 10 -5 , r 2 = 0.53) and OI (O:C = 0.35 + OI × 1.27 × 10 -3 , r 2 = 0.46). We used these equations to estimate the elemental H:C and O:C ratios of SOM, which allowed us to then calculate the SOC molar mass (MSOC, g SOC mol -1 SOM) using the following equation: = + : + : + : (4) where MC, MH, MO and MN are respectively the molar mass (g mol -1 ) of C, H, O and N.
Differential scanning calorimetry (DSC) during ramped combustion was also performed to measure the net energy released by SOM combustion (enthalpy of combustion), knowing that some of the energy applied to the sample is consumed by the breakdown of the organo-mineral associations. The mineral phase itself could also contribute to heat flux during ramped combustion, mainly generating endothermic reactions. However, previous studies showed that the contribution of the mineral phase to heat flux is usually small relatively to the contribution from SOM 28,29 . Briefly, ca. 50 mg of ground (<250 μm) samples were placed in 70 L alumina crucibles, with an identical empty crucible used as a reference, and subjected to oxidation ramping (25-1000 °C, ramping rate of 5 °C·min −1 under a synthetic CO2-free air atmosphere) using a thermal analyser simultaneously performing DSC and thermogravimetry analyses (TGA/DSC 3+ model, Mettler-Toledo, Greifensee, Switzerland). DSC heat fluxes (the exothermic or endothermic energy fluxes from the sample, referenced to an empty alumina crucible) were recorded every second, and the DSC thermograms were corrected a posteriori using a spline linear baseline. Net energy released was determined by integrating the exothermic region of the DSC thermogram (185-600 °C, Figure S3c), which represents the temperature range in which SOM is combusted 28 . Energy density of SOM (ΔE, in kJ g -1 SOC, Supplementary Table 8) was calculated as the net energy released divided by SOC content, which was recovered from thermogravimetry mass loss converted to SOC content with an equation accounting for mass loss due to clay water loss 30 . To get both ΔE and µEa expressed in the same unit, ΔE was converted in kJ mol -1 SOM by multiplying it with the SOC molar mass (g SOC mol -1 SOM) estimated based on C:H:O:N stoichiometry as explained above. We acknowledge that our SOC molar mass values remain approximative because of the uncertainty associated to the estimation of H:C and O:C ratios based on HI and OI indices. However, this allowed ΔE and µEa to be expressed in the same unit for the ROI calculation, and the estimated SOC molar mass varied little between treatments (Supplementary Table 8).

Isotopic partitioning
Correction of plant-soil system respiration for background atmospheric CO2 During the sampling of CO2 fluxes in the first series of incubations, the microcosms were sealed in opaque, airtight PVC chambers. Just before the sealing, each PVC chamber was intensively ventilated for 1 min with ambient air and we took care to avoid any contamination of the chamber air by breathing. The ambient air used for the ventilation was sampled to measure the initial amount and δ 13 C of CO2 in the chamber at the beginning of the incubation. Each microcosm sealed in chamber was then incubated for 24 h at temperature-controlled conditions (21.5 °C) in the laboratory. After 24 h of CO2 release by the plant-soil system, the chamber gas was sampled by transfer into a glass flask with a vacuum pump. After beeing ventilated ten times his volume the flask was airtight sealed until gas analysis. Its CO2 concentration as well as δ 13 C were measured using a Gas Chromatograph (Clarus 480, Perkin Elmer, Waltham, MA, USA) and an isotope laser spectrometer (G2201-i Isotopic Analyser, Picarro, Santa Clara, CA, USA). The amount and δ 13 C of CO2 derived from the plant-soil system respiration were corrected for background atmospheric CO2 using the following equations: where Rtotal and δ 13 Ctotal are the total amount and δ 13 C of CO2 release by the plant-soil system after correction; CO2-chamber and δ 13 Cchamber are the total amount and δ 13 C of CO2 measured in the chamber at the end of the incubation; and CO2-atmosphere and δ 13 Catmosphere are the total amount and δ 13 C of atmospheric CO2 sampled at the beginning of the incubation.

Plant-derived and soil-derived CO2 fluxes
In the first series of incubations, we observed that the CO2 derived from the respiration of the unplanted microcosms was systematically depleted in 13 C relative to the δ 13 C of SOC for most treatments (Supplementary Table 11), with an average differences in δ 13 C of -4.2 ‰, whereas CO2 derived from SOC respiration usually tend to be slightly 13 C enriched relative to SOC 31 . We concomitantly observed algae development on the soil surface of unplanted microcosms. Furthermore, Cros, et al. (2019) 32 observed the fixation of a small quantity of CO2 in unplanted soil linked to the development of algae during daytime. They also observed a depleted δ 13 C signature of unplanted soil respiration during 24 h dark incubations, which was of the same order of magnitude relative to their δ 13 C of SOC (~ -3 to -4 ‰). This highlights that a non-negligible photosynthetic activity occurred in the unplanted soil due to algae development, which could in turn respire labeled OC and explain the depleted δ 13 C of the unplanted soil respiration. Following the procedure used in a previous study 33 , we thus decided to apply the same isotopic partitioning on Rtotal of unplanted microcosms to correct for algae respiration, thus yielding more accurate estimation of Rsoil for the unplanted controls. The δ 13 C values of unplanted soil respiration measured at the end of experiment in the second series of incubations after removing the thin layer contaminated by algae on soil surface (~ 1 to 2 mm) were thus used as δ 13 Csoil for both planted and unplanted soil in the first series of incubations.
Cros et al.
(personal communication) also tested the possibility this depleted δ 13 C could be due to labelled CO2 back-diffusion during the 24 h dark incubation. After 24 h of pot ventilation either with labeled or with unlabeled (ambient) air, they measured the δ 13 C of the CO2 released during the incubation. In both cases, they found very similar δ 13 C, indicating that the contribution of labelled CO2 back-diffusion to the depleted δ 13 C signature of unplanted soil respiration was negligible. This is consistent with a recent study that also found very weak backdiffusion of labelled CO2 during incubation 34 .
Uncertainty in 13 C source partitioning related to sampling and analytical errors was calculated following Phillips & Gregg (2001) 35 with the following equation: where σfsoil is the standard error of the mean soil source proportion, that is fsoil; δ ̅ , δ ̅ , δ ̅ are respectively the mean δ 13 C of the mixture and of soil and plant sources; σ δ ̅ , σ δ ̅ and σ δ ̅ are respectively the standard error of mean δ 13 C of the mixture and of soil and plant sources; σ is the population standard deviation in δ 13 C among individual samples of source x; σ is the δ 13 C analytical standard deviation; and nx is the population size of source x. The σ values are respectively 0.10 and 0.12 ‰ for elemental analyser/isotope-ratio mass spectrometer and isotope laser spectrometer measurements. Since isotopic partitioning was performed at the microcosm level where the whole plant biomass has been sampled (plant source) and the isotopic signature has been measured on a well homogenized sample for both microcosm atmosphere (mixture) and plant biomass (plant source), the σ values for these two sources were assumed to be zero. We found that σfsoil values were on average 1.1 and 1.2 % respectively for the first and second series of incubations (Supplementary Tables 11 and 12), indicating low level of uncertainty associated with sampling and analytical errors 35 .
Our isotopic partitioning relied on several assumptions about the isotopic signature of our sources. For the first series of incubations, we used the mass-weighted δ 13 C of the mesocosm shoot and living root biomass as δ 13 Cplant, assuming negligible fractionation during whole-plant respiration 36 . We also used the mean δ 13 C of plant biomass across all treatments as δ 13 Cplant for unplanted soil, assuming similar δ 13 C fractionation during C3 photosynthesis for algae than for D. glomerata 37 . Finally, we used the δ 13 C values of unplanted control respiration measured at the end of experiment in the second series of incubations as δ 13 Csoil for both planted and unplanted soil in the first series of incubations, assuming constant δ 13 C fractionation of soilderived respiration across time. For the second series of incubations, we used δ 13 Cplant values taken as the living root biomass δ 13 C values corrected by the δ 13 C fractionation factor (f) of root respiration, assuming a constant f value of -0.61 ‰ 33 . We also used the Δ 14 C of root biomass as Δ 14 Cplant, assuming no fractionation of 14 C during respiration of root-derived OC.
In order to evaluate the uncertainty associated with our isotopic mixing model assumptions, we performed a sensitivity analysis where we quantified the error in fsoil and Δ 14 Csoil related to a 1 ‰ variation in δ 13 Csoil and δ 13 Cplant, and a 2 ‰ variation in Δ 14 Cplant. Given that the deviation can be positive or negative, we tested a deviation of +0.5 or +1 and −0.5 or −1 ‰ respectively corresponding to an amplitude of 1 or 2 ‰. The error was then expressed as the difference in fsoil and Δ 14 Csoil between +0.5 or +1 and −0.5 or −1 ‰ source values. We found that the average errors in fsoil related to a 1 ‰ variation in δ 13 Cplant was respectively 2.0 and 1.9 % for the first and second series of incubations, while the average error related to a 1 ‰ variation in δ 13 Csoil was 1.5 % for the first series of incubations (Supplementary Tables 11 and 12). We also found that the average errors in the absolute values of RPE related to a 1 ‰ variation in δ 13 Cplant was respectively 19.3 and 12.2 % for the first and second series of incubations, while the average error related to a 1 ‰ variation in δ 13 Csoil was 0.2 % for the first series of incubations (Supplementary Tables 11 and 12). The average errors in Δ 14 Csoil related to a 1 ‰ variation in δ 13 Cplant and a 2 ‰ variation Δ 14 Cplant were respectively 19.4 and 2.8 ‰ (Supplementary Table  13). These low levels of uncertainty provide evidence that our isotopic partitioning was robust.
For the first and second series of incubations, the average δ 13 C difference between plant and soil sources were respectively 24.9 and 24.8 ‰ (range of respectively 23.6 to 26.4, and 22.9 to 26.6 across treatments). This is substantially larger than source difference typical yielded by the natural 13 C-labeling method 38 , as well as the traditional continuous 13 C-labeling method based on adding fossil fuel-derived CO2 in a flow of ambient air 39 . This strong labeling allows to improve the accuracy of isotopic partitioning 31,32,35 .

Living root biomass
As the root material harvested for topsoil after the second incubation series was composed of both pre-existing root litter (unlabelled) and living root (labelled) that could not be clearly visually sorted, we used an isotopic partitioning method to estimate the biomass of living roots for each planted topsoil core. Assuming negligible pre-existing (unlabeled) root litter biomass in subsoil at the end of the experiment and equal δ 13 C difference between shoot and root biomass (Δδ 13 Cshoot-root) in topsoil and subsoil, we calculated the δ 13 C of living root biomass as the microcosm shoot δ 13 C minus the Δδ 13 Cshoot-root of subsoil (Supplementary Table 10). Further assuming equal δ 13 C for pre-existing (dead) root litter and SOC, we used the following equation based on a two-source isotopic mixing model.
where Roottotal and δ 13 Ctotal are respectively the biomass and δ 13 C of both dead and living roots; Rootliving and δ 13 Cliving are respectively the biomass and δ 13 C of both living roots; and δ 13 Cdead is the δ 13 C of dead roots. The standard errors in 13 C source partitioning related to sampling and analytical errors was on average 0.6 %, indicating low uncertainty level (Supplementary Table 10). A similar sensitivity analysis performed as described above showed that the average errors in 13 C source partitioning related to a 1 ‰ variation in δ 13 Cliving was 3.7 %, while the average error related to a 1 ‰ variation in δ 13 Cdead was 0.6 % (Supplementary Table 10). These low levels of uncertainty provide evidence that our isotopic partitioning was here also robust.

Net rhizodeposition
We quantified the net rhizodeposition corresponding to the root-derived SOC remaining in the soil after microbial utilization. Net rhizodeposition was used here as a proxy of the gross rhizodeposition corresponding to the flux of fresh OC supply into the soil via living roots, which remains so far very challenging to quantify 40 . We acknowledge that net rhizodeposition is not only driven by gross rhizodeposition, but also by other factors affecting the stabilization and destabilization of rhizodeposits such as soil mineralogy or microbial communities. However, we argue that net rhizodeposition remains a good proxy of gross rhizodeposition at the microcosm scale, and provides here useful insights about how SOC decomposition response to variation in living root density with core depth in the second incubation series was related to variation in gross rhizodeposition.
Net rhizodeposition was estimated for each planted soil core harvested after the second incubation series using the following equation based on a two-source isotopic mixing model: where SOCtotal and δ 13 CSOC-final are respectively the SOC content and δ 13 C of the planted soil core at the end of the experiment, δ 13 CSOC-initial is the average δ 13 C of SOC from soil cores sampled in the field at the beginning of the experiment, and δ 13 Cplant is the δ 13 C of living root biomass of the planted microcosm. A sensitivity analysis was performed as described above to assess the uncertainty associated to our assumption of a negligible 13 C fractionation of root-derived OC during microbial utilization. The average errors in 13 C source partitioning related to a 1 ‰ variation in δ 13 Croot was only 0.06 %, indicating low levels of uncertainty related to this assumption (Supplementary  Table 14). However, the standard error in 13 C source partitioning related to sampling and analytical errors was on average 0.77 %, which we acknowledge is rather high relative the average proportion of root-derived SOC found at the end of the experiment (1.37 %, Supplementary Table 14). Such level of uncertainty nevertheless remains within the range of previous experiments quantifying plant-derived OC in a large reservoir of existing SOC with a comparable labeling intensity 33,[41][42][43] .

Statistical analyses
We used a rotated principal component analysis to explore soil properties covariance and divergence divergence between treatments (Fig. 1a). The rotated principal component analysis (RCA) was performed using the 'principal' function of the 'Psych' package 44 . Rotation is commonly used in principal component analysis to simplify interpretation of principal components by maximizing/minimizing the correlations between factors and component axes. In order to simplify the RCA ordination, we selected only two axes. Analyses of variance were used to partition the variance explained by the factors 'soil layer', 'soil type' and their interaction in the two first axis scores (Fig. 1b) and soil properties (Supplementary Table 2).
Partial η² of depth effect on SOC properties is calculated as the sum of squares for the depth effect divided by the total sum of squares (after accounting for the variance associated with soil type effect). It was computed using the 'eta_squared' function of the 'effectsize' package 45 on a linear mixed-effect model fitted using the 'lmer' function of the 'lme4' package 46 and including 'soil layer' as a fixed factor and 'soil type' as a random factor.
For each series of incubations, the responses of kSOC, RPE and Δ 14 Csoil to predictors were assessed by ordinary least squares regression for each treatment (Figs. 3, 4 and 5). We tested linear (Y = a + bX), polynomial (Y = a + bX + cX 2 ) and power (Y = aX b ) regression functions, where Y is the response variable and X is the predictor. The kSOC and RPE values were standardized for each treatment to a common high value of the following predictors: 'respiration of plant-derived OC' for the first incubation series, as well as 'living root density', 'respiration of plant-derived OC' and 'net rhizodeposition' for the second incubation series. To do so, we used the regression model parameters to predict kSOC and RPE values at the mean predictor value across treatments of the last sampling time in the first incubation series, corresponding to 9.94 g C-CO2 m -2 day -1 for 'respiration of plant-derived OC'. For the second incubation series, the same procedure was applied with the mean predictor values across treatments in the 0-20 cm depth soil core, corresponding to 3.57 g dm -3 for 'living root density', 20.48 mg C-CO2 dm -3 day -1 for 'respiration of plant-derived OC' and 0.79 g SOCroot dm -3 for 'net rhizodeposition'. Additionally, we used analyses of covariance including the quantitative explanatory variables, the factors 'soil layers' and 'soil type', and their interactions as fixed factors to test their effects and quantify the proportion of variance they explain (Supplementary Tables 3, 4 and 5). To deal with the repeated measures design in both the first and second incubation series, we used linear mixed-effect models including 'microcosm' as random factor in regression and analyses of covariance. Linear mixed-effects models were fitted using the 'lmer' function of the 'lme4' package 46 . Statistical significance of fixed predictors were assessed based on Satterthwaite's approximation of denominator degrees of freedom using the 'anova' function of the 'lmerTest' package 47 . Marginal and conditional r 2 were computed based on Nakagawa and Schielzeth 48 using the 'r.squaredGLMM' function of the 'MuMIn' package 49 .
To examine the relationship between soil biogeochemistry and organic matter dynamics, we explore partial bivariate relationships between variation in SOC dynamics and SOM properties across depth while controlling for soil types. Before quantification of bivariate relationships, each variable was centered using the residuals of a linear model with 'soil type' as a fixed effect. Partial bivariate relationships of radiocarbon-based mean SOC age with both the return-onenergy-investment of microbial SOC decomposition (ROI) and root density were first examined using ordinary least squares linear regression models (Fig. 2). The slope coefficients of the regression models were standardized by range. To do so, the unstandardized slope coefficient was multiped by the range (the difference between the maximum and minimum values) of the predictor and divided by the range of the dependant variable). It was computed using the 'coefs' function of the 'piecewiseSEM' package 50 Partial correlations between SOC dynamics variables, including 14 C-based SOC turnover time as well as unplanted kSOC and standardized RPE of each incubation series, and SOM properties were then also performed by computing Spearman's correlation coefficients ( Supplementary Fig. 3). It was computed using the 'rcorr' function of the 'Hmisc' package 51 with depth-centered variables. Additionally, we performed an ordination of the same SOC dynamics variables constrained by soil biogeochemical predictors (same set of variables used in the RCA) using a redundancy analysis ( Supplementary  Fig. 4). The redundancy analysis (RDA) was performed using the 'rda' function of the 'vegan' package 52 with raw variables. The significance of the overall constrained ordination and of the two first axes were tested by permutation tests 53 , using the 'anova' function of the 'vegan' package 52 (1,000 permutations).
All analyses were performed using R v3.4.3 54 . Null hypothesis testing was always based on two-sided statistical tests. The normal distribution and homogeneity of variances of the model residuals were graphically checked and data were log-transformed when necessary.

Supplementary Figures
Supplementary Fig. 1. Energetic properties of soil organic matter across treatments. a, 3 thermogram of soil organic carbon (SOC) decomposition measured by Rock-Eval sequential 4 pyrolysis and oxidation. b, distribution of activation energy, p(Ea). c, thermogram of heat flow 5 measured by differential scanning calorimetry. Polygons around lines represent 95% 6 confidence intervals around treatment means (n = 4 replicate soil cores). 7 10 8 Supplementary Fig. 2. Soil core sampling and microcosm preparation method. A percussion core drill equipped with a steel tube that can be 9 opened from sideways was used to extract intact soil columns of 8 cm diameter (a, b, d). For each layer, three soil cores collected with a knife (c) 10 from the same depth in the initial soil column were gently stacked together (a, d), tightly sealed within a polyethylene sheath (d, e) and transferred 11 into a bottom-capped polyvinyl chloride (PVC) tube (e) to form a new soil column exclusively made of topsoil or subsoil (a). See Supplementary 12 Fig. 5 for the soil core sampling design of each soil type. 13 11

14
Supplementary Fig. 3. Heatmap of partial correlations between variation in soil organic carbon (SOC) dynamics parameters and soil organic 15 matter (SOM) properties across depth (n = 24 soil cores/microcosms). Each variable was normalized for variation across soil type using the residuals 16 of a linear model with 'soil type' as a fixed effect. The variation in native SOC decomposition rate (kSOC) of unplanted soil and rhizosphere priming 17 effect (RPE) across time and soil column depth respectively represent the parameters from the first and second series of incubations. RPE values 18 were standardized to a common high predictor value across treatments. Root density values from the biogeochemical characterization of soil are 19 used here to reflect in situ root density. See Supplementary Fig. 1 for abbreviations. ***, P < 0.001; **, P < 0.01; *, P < 0.05; †, P < 0. 10 Subsoil Topsoil Subsoil Topsoil Subsoil SOC molar mass (g SOC mol -1 SOM) 6.11 ± 0.02 5.86 ± 0.01 6.51 ± 0.01 6.59 ± 0.01 6.28 ± 0.00 6.09 ± 0.06 Energy density (kJ g -1 SOC) 25 13.0 ± 1.7 13.8 ± 0.3 11.9 ± 0.6 11.9 ± 0.2 19.9 ± 3.0 25.2 ± 0.5 Living plant N content (mg kg -1 ) 6.0 ± 0.4 5.4 ± 0.3 6.7 ± 0.3 5.0 ± 0.5 6.7 ± 0.3 7.4 ± 0.8 ⁋ Post-fertilization. fliving is the proportion of living roots in the living and dead root mixture; σfliving is the standard error of fliving related to sampling and analytical errors; Δfliving|Δδ 13 Cliving and Δfliving|Δδ 13 Cdead are the variations in fliving to a 1 ‰ variation in δ 13 Cliving and δ 13 Cdead, respectively.  13 Csoil and δ 13 Cplant are the δ 13 C of respectively total, native soil organic carbon and root-derived organic carbon respiration. fsoil is the proportion of CO2 from native SOC respiration; σfsoil is the standard error of fsoil related to sampling and analytical errors; Δfsoil|Δδ 13 Cplant and ΔRPE|Δδ 13 Cplant are the variation in respectively fsoil and RPE to a 1 ‰ variation in δ 13 Cplant.  13 Ctotal, δ 13 Csoil and δ 13 Cplant are the δ 13 C of respectively total, native and root-derived soil organic carbon. froot is the proportion of root-derived soil organic carbon; σfroot is the standard error of froot related to sampling and analytical errors; Δfroot|Δδ 13 Croot is the variation in fsoil to a 1 ‰ variation in δ 13 Croot.